***Set up and run simulations***
********************************
//run "prepare_datasets.do" to recreate all datasets used in analyses from original data
clear all
cap program drop psim pdata
set rmsg on
set seed 5000
//install programs if needed
//ssc install reghdfe,replace
//ssc install ftools,replace 
//ssc install regsave,replace
//ssc install tsspell,replace

*Set up programs for simulating panel datasets*
//dataset simulation for main results, restricted time period and number of countries, and specifications that include lagged dependent variable
program define pdata
    drop _all
    set obs 225
    gen id = _n
foreach i in 2 3 4 5 6 7 8 9 10 11 12 13 14 {
	//randomly combine time-series of democracy data and GDP per capita data for different countries
g temp=runiform()
	sort temp
	g id`i'=_n
	drop temp	
}
//drop observations that match actual countries for both time-series
drop if id==id2
drop if id3==id4
drop if id5==id6
drop if id7==id8
drop if id9==id10
drop if id11==id12
drop if id13==id14
	
    expand 250
	save "base.dta",replace
	g d=1
	bys id: gen t=sum(d)
	drop d
	replace t=t+1788
qui	merge 1:1 id t using "democracy_for_simu.dta"
	keep if _merge==3
	drop _merge
qui	merge 1:1 id2 t using "gdp_for_simu.dta"
	keep if _merge==3
	drop _merge	
	save "temp.dta",replace
	use "base.dta",replace
	g d=1
	bys id3: gen t=sum(d)
	drop d
	replace t=t+1788
	drop id
	rename id3 id
qui	merge 1:1 id t using "democracy_for_simu.dta"
	keep if _merge==3
	
	drop _merge
	drop id2
	rename id4 id2
qui	merge 1:1 id2 t using "gdp_for_simu.dta"
	keep if _merge==3
	
	drop _merge
	replace id=id+300
	save "temp2.dta",replace
	use "base.dta",replace
	g d=1
	bys id5: gen t=sum(d)
	drop d
	replace t=t+1788
	drop id
	rename id5 id
qui	merge 1:1 id t using "democracy_for_simu.dta"
	keep if _merge==3

	drop _merge
	drop id2
	rename id6 id2
qui	merge 1:1 id2 t using "gdp_for_simu.dta"
	keep if _merge==3

	drop _merge
	replace id=id+600
	save "temp3.dta",replace
	use "base.dta",replace
	g d=1
	bys id5: gen t=sum(d)
	drop d
	replace t=t+1788
	drop id
	rename id7 id
qui	merge 1:1 id t using "democracy_for_simu.dta"
	keep if _merge==3

	drop _merge
	drop id2
	rename id8 id2
qui	merge 1:1 id2 t using "gdp_for_simu.dta"
	keep if _merge==3

	drop _merge
	replace id=id+900
	save "temp4.dta",replace
	use "base.dta",replace
	g d=1
	bys id5: gen t=sum(d)
	drop d
	replace t=t+1788
	drop id
	rename id9 id
qui	merge 1:1 id t using "democracy_for_simu.dta"
	keep if _merge==3

	drop _merge
	drop id2
	rename id10 id2
qui	merge 1:1 id2 t using "gdp_for_simu.dta"
	keep if _merge==3

	drop _merge
	replace id=id+1200
	save "temp5.dta",replace
		use "base.dta",replace
	g d=1
	bys id5: gen t=sum(d)
	drop d
	replace t=t+1788
	drop id
	rename id11 id
qui	merge 1:1 id t using "democracy_for_simu.dta"
	keep if _merge==3

	drop _merge
	drop id2
	rename id12 id2
qui	merge 1:1 id2 t using "gdp_for_simu.dta"
	keep if _merge==3

	drop _merge
	replace id=id+1500
	save "temp6.dta",replace
		use "base.dta",replace
	g d=1
	bys id5: gen t=sum(d)
	drop d
	replace t=t+1788
	drop id
	rename id13 id
qui	merge 1:1 id t using "democracy_for_simu.dta"
	keep if _merge==3

	drop _merge
	drop id2
	rename id14 id2
qui	merge 1:1 id2 t using "gdp_for_simu.dta"
	keep if _merge==3

	drop _merge
	replace id=id+1800
	save "temp7.dta",replace
append using "temp2.dta"
	append using "temp.dta"
	append using "temp3.dta"
	append using "temp4.dta"
	append using "temp5.dta"
	append using "temp6.dta"
	drop if id==.
	save "tempfinal.dta",replace
	duplicates drop id,force
	g rid=rnormal(0,100)
	egen rid2=group(rid)
	keep id rid2
	merge 1:m id using "tempfinal.dta"
	drop _merge id
	rename rid2 id
	sort id t
	//Ensure that simulated dataset have a similar number of countries as observed in the population for different time periods
	g keep=0
	g tid=rnormal(0,100) if t==1800
	egen k=group(tid) if t==1800
	by id: egen maxk=max(k)
	recode keep (0=1) if maxk<16 & maxk!=. & t>1799
	drop k maxk tid
	
	g tid=rnormal(0,100) if t==1816
	egen k=group(tid) if t==1816 & keep!=1
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1816
	egen obs=sum(d)
recode keep (0=1) if maxk<(31-obs) & maxk!=. & t>1815
	drop k maxk tid obs d
	g tid=rnormal(0,100) if t==1820
	egen k=group(tid) if t==1820 & keep!=1
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1820
	egen obs=sum(d)
	recode keep (0=1) if maxk<(38-obs) & maxk!=. & t>1819
	drop k maxk tid obs d
	g tid=rnormal(0,100) if t==1825
	egen k=group(tid) if t==1825 & keep!=1
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1825
	egen obs=sum(d)
	recode keep (0=1) if maxk<(42-obs) & maxk!=. & t>1824
	drop k maxk tid obs d
g tid=rnormal(0,100) if t==1830
	egen k=group(tid) if t==1830 & keep!=1
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1830
	egen obs=sum(d)
	recode keep (0=1) if maxk<(46-obs) & maxk!=. & t>1829
	drop k maxk tid obs d
	g tid=rnormal(0,100) if t==1838 & keep!=1
	egen k=group(tid) if t==1838
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1838
	egen obs=sum(d)
	recode keep (0=1) if maxk<(49-obs) & maxk!=. & t>1837
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1844 & keep!=1
	egen k=group(tid) if t==1844
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1844
	egen obs=sum(d)
	recode keep (0=1) if maxk<(53-obs) & maxk!=. & t>1843
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1857 & keep!=1
	egen k=group(tid) if t==1857
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1857
	egen obs=sum(d)
	recode keep (0=1) if maxk<(54-obs) & maxk!=. & t>1856
	drop k maxk tid obs d
	g tid=rnormal(0,100) if t==1901 & keep!=1
	egen k=group(tid) if t==1901
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1901
	egen obs=sum(d)
	recode keep (0=1) if maxk<(57-obs) & maxk!=. & t>1900
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1918 & keep!=1
	egen k=group(tid) if t==1918
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1918
	egen obs=sum(d)
	recode keep (0=1) if maxk<(67-obs) & maxk!=. & t>1917
	drop k maxk tid obs d
g tid=rnormal(0,100) if t==1945 & keep!=1
	egen k=group(tid) if t==1945
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1945
	egen obs=sum(d)
	recode keep (0=1) if maxk<(71-obs) & maxk!=. & t>1944
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1948 & keep!=1
	egen k=group(tid) if t==1948
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1948
	egen obs=sum(d)
	recode keep (0=1) if maxk<(82-obs) & maxk!=. & t>1947
	drop k maxk tid obs d
	g tid=rnormal(0,100) if t==1956 & keep!=1
	egen k=group(tid) if t==1956
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1956
	egen obs=sum(d)
	recode keep (0=1) if maxk<(91-obs) & maxk!=. & t>1955
	drop k maxk tid obs d
	g tid=rnormal(0,100) if t==1960 & keep!=1
	egen k=group(tid) if t==1960
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1960
	egen obs=sum(d)
	recode keep (0=1) if maxk<(111-obs) & maxk!=. & t>1959
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1962 & keep!=1
	egen k=group(tid) if t==1962
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1962
	egen obs=sum(d)
	recode keep (0=1) if maxk<(121-obs) & maxk!=. & t>1961
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1966 & keep!=1
	egen k=group(tid) if t==1966
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1966
	egen obs=sum(d)
	recode keep (0=1) if maxk<(136-obs) & maxk!=. & t>1965
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1971 & keep!=1
	egen k=group(tid) if t==1971
	by id: egen maxk=max(k) 
	g d=1 if keep==1 & t==1971
	egen obs=sum(d)
	recode keep (0=1) if maxk<(141-obs) & maxk!=. & t>1970
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1975 & keep!=1
	egen k=group(tid) if t==1975
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1975
	egen obs=sum(d)
	recode keep (0=1) if maxk<(151-obs) & maxk!=. & t>1974
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1991 & keep!=1
	egen k=group(tid) if t==1991
	
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1991
	egen obs=sum(d)
	recode keep (0=1) if maxk<(171-obs) & maxk!=. & t>1990
	drop if keep!=1
	keep id t x lng time_* 
	save "tempfinal.dta",replace
end
//dataset simulation with civil war onset as outcome
program define pdata2
    drop _all
    set obs 225
    gen id = _n
foreach i in 2 3 4 5 6 7 8 9 10 11 12 13 14 {
g temp=runiform()
	sort temp
	g id`i'=_n
	drop temp	
}
drop if id==id2
drop if id3==id4
drop if id5==id6
drop if id7==id8
drop if id9==id10
drop if id11==id12
drop if id13==id14
	
    expand 250
	save "base.dta",replace
	g d=1
	bys id: gen t=sum(d)
	drop d
	replace t=t+1788
qui	merge 1:1 id t using "democracy_for_simu.dta"
	keep if _merge==3
	drop _merge
qui	merge 1:1 id2 t using "cw_onset.dta"
	keep if _merge==3
	drop _merge	
	save "temp.dta",replace
	use "base.dta",replace
	g d=1
	bys id3: gen t=sum(d)
	drop d
	replace t=t+1788
	drop id
	rename id3 id
qui	merge 1:1 id t using "democracy_for_simu.dta"
	keep if _merge==3
	
	drop _merge
	drop id2
	rename id4 id2
qui	merge 1:1 id2 t using "cw_onset.dta"
	keep if _merge==3
	
	drop _merge
	replace id=id+300
	save "temp2.dta",replace
	use "base.dta",replace
	g d=1
	bys id5: gen t=sum(d)
	drop d
	replace t=t+1788
	drop id
	rename id5 id
qui	merge 1:1 id t using "democracy_for_simu.dta"
	keep if _merge==3

	drop _merge
	drop id2
	rename id6 id2
qui	merge 1:1 id2 t using "cw_onset.dta"
	keep if _merge==3

	drop _merge
	replace id=id+600
	save "temp3.dta",replace
	use "base.dta",replace
	g d=1
	bys id5: gen t=sum(d)
	drop d
	replace t=t+1788
	drop id
	rename id7 id
qui	merge 1:1 id t using "democracy_for_simu.dta"
	keep if _merge==3

	drop _merge
	drop id2
	rename id8 id2
qui	merge 1:1 id2 t using "cw_onset.dta"
	keep if _merge==3

	drop _merge
	replace id=id+900
	save "temp4.dta",replace
	use "base.dta",replace
	g d=1
	bys id5: gen t=sum(d)
	drop d
	replace t=t+1788
	drop id
	rename id9 id
qui	merge 1:1 id t using "democracy_for_simu.dta"
	keep if _merge==3

	drop _merge
	drop id2
	rename id10 id2
qui	merge 1:1 id2 t using "cw_onset.dta"
	keep if _merge==3

	drop _merge
	replace id=id+1200
	save "temp5.dta",replace
		use "base.dta",replace
	g d=1
	bys id5: gen t=sum(d)
	drop d
	replace t=t+1788
	drop id
	rename id11 id
qui	merge 1:1 id t using "democracy_for_simu.dta"
	keep if _merge==3

	drop _merge
	drop id2
	rename id12 id2
qui	merge 1:1 id2 t using "cw_onset.dta"
	keep if _merge==3

	drop _merge
	replace id=id+1500
	save "temp6.dta",replace
		use "base.dta",replace
	g d=1
	bys id5: gen t=sum(d)
	drop d
	replace t=t+1788
	drop id
	rename id13 id
qui	merge 1:1 id t using "democracy_for_simu.dta"
	keep if _merge==3

	drop _merge
	drop id2
	rename id14 id2
qui	merge 1:1 id2 t using "cw_onset.dta"
	keep if _merge==3

	drop _merge
	replace id=id+1800
	save "temp7.dta",replace
append using "temp2.dta"
	append using "temp.dta"
	append using "temp3.dta"
	append using "temp4.dta"
	append using "temp5.dta"
	append using "temp6.dta"
	drop if id==.
	save "tempfinal.dta",replace
	duplicates drop id,force
	g rid=rnormal(0,100)
	egen rid2=group(rid)
	keep id rid2
	merge 1:m id using "tempfinal.dta"
	drop _merge id
	rename rid2 id
	sort id t
	g keep=0
	g tid=rnormal(0,100) if t==1800
	egen k=group(tid) if t==1800
	by id: egen maxk=max(k)
	recode keep (0=1) if maxk<16 & maxk!=. & t>1799
	drop k maxk tid
	
	g tid=rnormal(0,100) if t==1816
	egen k=group(tid) if t==1816 & keep!=1
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1816
	egen obs=sum(d)
recode keep (0=1) if maxk<(31-obs) & maxk!=. & t>1815
	drop k maxk tid obs d
	g tid=rnormal(0,100) if t==1820
	egen k=group(tid) if t==1820 & keep!=1
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1820
	egen obs=sum(d)
	recode keep (0=1) if maxk<(38-obs) & maxk!=. & t>1819
	drop k maxk tid obs d
	g tid=rnormal(0,100) if t==1825
	egen k=group(tid) if t==1825 & keep!=1
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1825
	egen obs=sum(d)
	recode keep (0=1) if maxk<(42-obs) & maxk!=. & t>1824
	drop k maxk tid obs d
g tid=rnormal(0,100) if t==1830
	egen k=group(tid) if t==1830 & keep!=1
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1830
	egen obs=sum(d)
	recode keep (0=1) if maxk<(46-obs) & maxk!=. & t>1829
	drop k maxk tid obs d
	g tid=rnormal(0,100) if t==1838 & keep!=1
	egen k=group(tid) if t==1838
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1838
	egen obs=sum(d)
	recode keep (0=1) if maxk<(49-obs) & maxk!=. & t>1837
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1844 & keep!=1
	egen k=group(tid) if t==1844
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1844
	egen obs=sum(d)
	recode keep (0=1) if maxk<(53-obs) & maxk!=. & t>1843
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1857 & keep!=1
	egen k=group(tid) if t==1857
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1857
	egen obs=sum(d)
	recode keep (0=1) if maxk<(54-obs) & maxk!=. & t>1856
	drop k maxk tid obs d
	g tid=rnormal(0,100) if t==1901 & keep!=1
	egen k=group(tid) if t==1901
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1901
	egen obs=sum(d)
	recode keep (0=1) if maxk<(57-obs) & maxk!=. & t>1900
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1918 & keep!=1
	egen k=group(tid) if t==1918
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1918
	egen obs=sum(d)
	recode keep (0=1) if maxk<(67-obs) & maxk!=. & t>1917
	drop k maxk tid obs d
g tid=rnormal(0,100) if t==1945 & keep!=1
	egen k=group(tid) if t==1945
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1945
	egen obs=sum(d)
	recode keep (0=1) if maxk<(71-obs) & maxk!=. & t>1944
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1948 & keep!=1
	egen k=group(tid) if t==1948
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1948
	egen obs=sum(d)
	recode keep (0=1) if maxk<(82-obs) & maxk!=. & t>1947
	drop k maxk tid obs d
	g tid=rnormal(0,100) if t==1956 & keep!=1
	egen k=group(tid) if t==1956
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1956
	egen obs=sum(d)
	recode keep (0=1) if maxk<(91-obs) & maxk!=. & t>1955
	drop k maxk tid obs d
	g tid=rnormal(0,100) if t==1960 & keep!=1
	egen k=group(tid) if t==1960
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1960
	egen obs=sum(d)
	recode keep (0=1) if maxk<(111-obs) & maxk!=. & t>1959
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1962 & keep!=1
	egen k=group(tid) if t==1962
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1962
	egen obs=sum(d)
	recode keep (0=1) if maxk<(121-obs) & maxk!=. & t>1961
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1966 & keep!=1
	egen k=group(tid) if t==1966
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1966
	egen obs=sum(d)
	recode keep (0=1) if maxk<(136-obs) & maxk!=. & t>1965
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1971 & keep!=1
	egen k=group(tid) if t==1971
	by id: egen maxk=max(k) 
	g d=1 if keep==1 & t==1971
	egen obs=sum(d)
	recode keep (0=1) if maxk<(141-obs) & maxk!=. & t>1970
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1975 & keep!=1
	egen k=group(tid) if t==1975
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1975
	egen obs=sum(d)
	recode keep (0=1) if maxk<(151-obs) & maxk!=. & t>1974
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1991 & keep!=1
	egen k=group(tid) if t==1991
	
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1991
	egen obs=sum(d)
	recode keep (0=1) if maxk<(171-obs) & maxk!=. & t>1990
	drop if keep!=1
	keep id t x onset
	save "tempfinal.dta",replace
end	
//dataset simulation with polyarchy index as treatment
program define pdata3
    drop _all
    set obs 225
    gen id = _n
foreach i in 2 3 4 5 6 7 8 9 10 11 12 13 14 {
g temp=runiform()
	sort temp
	g id`i'=_n
	drop temp	
}
drop if id==id2
drop if id3==id4
drop if id5==id6
drop if id7==id8
drop if id9==id10
drop if id11==id12
drop if id13==id14
	
    expand 250
	save "base.dta",replace
	g d=1
	bys id: gen t=sum(d)
	drop d
	replace t=t+1788
qui	merge 1:1 id t using "democracy_for_simu_poly.dta"
	keep if _merge==3
	drop _merge
qui	merge 1:1 id2 t using "gdp_for_simu.dta"
	keep if _merge==3
	drop _merge	
	save "temp.dta",replace
	use "base.dta",replace
	g d=1
	bys id3: gen t=sum(d)
	drop d
	replace t=t+1788
	drop id
	rename id3 id
qui	merge 1:1 id t using "democracy_for_simu_poly.dta"
	keep if _merge==3
	
	drop _merge
	drop id2
	rename id4 id2
qui	merge 1:1 id2 t using "gdp_for_simu.dta"
	keep if _merge==3
	
	drop _merge
	replace id=id+300
	save "temp2.dta",replace
	use "base.dta",replace
	g d=1
	bys id5: gen t=sum(d)
	drop d
	replace t=t+1788
	drop id
	rename id5 id
qui	merge 1:1 id t using "democracy_for_simu_poly.dta"
	keep if _merge==3

	drop _merge
	drop id2
	rename id6 id2
qui	merge 1:1 id2 t using "gdp_for_simu.dta"
	keep if _merge==3

	drop _merge
	replace id=id+600
	save "temp3.dta",replace
	use "base.dta",replace
	g d=1
	bys id5: gen t=sum(d)
	drop d
	replace t=t+1788
	drop id
	rename id7 id
qui	merge 1:1 id t using "democracy_for_simu_poly.dta"
	keep if _merge==3

	drop _merge
	drop id2
	rename id8 id2
qui	merge 1:1 id2 t using "gdp_for_simu.dta"
	keep if _merge==3

	drop _merge
	replace id=id+900
	save "temp4.dta",replace
	use "base.dta",replace
	g d=1
	bys id5: gen t=sum(d)
	drop d
	replace t=t+1788
	drop id
	rename id9 id
qui	merge 1:1 id t using "democracy_for_simu_poly.dta"
	keep if _merge==3

	drop _merge
	drop id2
	rename id10 id2
qui	merge 1:1 id2 t using "gdp_for_simu.dta"
	keep if _merge==3

	drop _merge
	replace id=id+1200
	save "temp5.dta",replace
		use "base.dta",replace
	g d=1
	bys id5: gen t=sum(d)
	drop d
	replace t=t+1788
	drop id
	rename id11 id
qui	merge 1:1 id t using "democracy_for_simu_poly.dta"
	keep if _merge==3

	drop _merge
	drop id2
	rename id12 id2
qui	merge 1:1 id2 t using "gdp_for_simu.dta"
	keep if _merge==3

	drop _merge
	replace id=id+1500
	save "temp6.dta",replace
		use "base.dta",replace
	g d=1
	bys id5: gen t=sum(d)
	drop d
	replace t=t+1788
	drop id
	rename id13 id
qui	merge 1:1 id t using "democracy_for_simu_poly.dta"
	keep if _merge==3

	drop _merge
	drop id2
	rename id14 id2
qui	merge 1:1 id2 t using "gdp_for_simu.dta"
	keep if _merge==3

	drop _merge
	replace id=id+1800
	save "temp7.dta",replace
append using "temp2.dta"
	append using "temp.dta"
	append using "temp3.dta"
	append using "temp4.dta"
	append using "temp5.dta"
	append using "temp6.dta"
	drop if id==.
	save "tempfinal.dta",replace
	duplicates drop id,force
	g rid=rnormal(0,100)
	egen rid2=group(rid)
	keep id rid2
	merge 1:m id using "tempfinal.dta"
	drop _merge id
	rename rid2 id
	sort id t
	g keep=0
	g tid=rnormal(0,100) if t==1800
	egen k=group(tid) if t==1800
	by id: egen maxk=max(k)
	recode keep (0=1) if maxk<16 & maxk!=. & t>1799
	drop k maxk tid
	
	g tid=rnormal(0,100) if t==1816
	egen k=group(tid) if t==1816 & keep!=1
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1816
	egen obs=sum(d)
recode keep (0=1) if maxk<(31-obs) & maxk!=. & t>1815
	drop k maxk tid obs d
	g tid=rnormal(0,100) if t==1820
	egen k=group(tid) if t==1820 & keep!=1
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1820
	egen obs=sum(d)
	recode keep (0=1) if maxk<(38-obs) & maxk!=. & t>1819
	drop k maxk tid obs d
	g tid=rnormal(0,100) if t==1825
	egen k=group(tid) if t==1825 & keep!=1
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1825
	egen obs=sum(d)
	recode keep (0=1) if maxk<(42-obs) & maxk!=. & t>1824
	drop k maxk tid obs d
g tid=rnormal(0,100) if t==1830
	egen k=group(tid) if t==1830 & keep!=1
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1830
	egen obs=sum(d)
	recode keep (0=1) if maxk<(46-obs) & maxk!=. & t>1829
	drop k maxk tid obs d
	g tid=rnormal(0,100) if t==1838 & keep!=1
	egen k=group(tid) if t==1838
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1838
	egen obs=sum(d)
	recode keep (0=1) if maxk<(49-obs) & maxk!=. & t>1837
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1844 & keep!=1
	egen k=group(tid) if t==1844
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1844
	egen obs=sum(d)
	recode keep (0=1) if maxk<(53-obs) & maxk!=. & t>1843
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1857 & keep!=1
	egen k=group(tid) if t==1857
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1857
	egen obs=sum(d)
	recode keep (0=1) if maxk<(54-obs) & maxk!=. & t>1856
	drop k maxk tid obs d
	g tid=rnormal(0,100) if t==1901 & keep!=1
	egen k=group(tid) if t==1901
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1901
	egen obs=sum(d)
	recode keep (0=1) if maxk<(57-obs) & maxk!=. & t>1900
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1918 & keep!=1
	egen k=group(tid) if t==1918
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1918
	egen obs=sum(d)
	recode keep (0=1) if maxk<(67-obs) & maxk!=. & t>1917
	drop k maxk tid obs d
g tid=rnormal(0,100) if t==1945 & keep!=1
	egen k=group(tid) if t==1945
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1945
	egen obs=sum(d)
	recode keep (0=1) if maxk<(71-obs) & maxk!=. & t>1944
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1948 & keep!=1
	egen k=group(tid) if t==1948
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1948
	egen obs=sum(d)
	recode keep (0=1) if maxk<(82-obs) & maxk!=. & t>1947
	drop k maxk tid obs d
	g tid=rnormal(0,100) if t==1956 & keep!=1
	egen k=group(tid) if t==1956
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1956
	egen obs=sum(d)
	recode keep (0=1) if maxk<(91-obs) & maxk!=. & t>1955
	drop k maxk tid obs d
	g tid=rnormal(0,100) if t==1960 & keep!=1
	egen k=group(tid) if t==1960
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1960
	egen obs=sum(d)
	recode keep (0=1) if maxk<(111-obs) & maxk!=. & t>1959
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1962 & keep!=1
	egen k=group(tid) if t==1962
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1962
	egen obs=sum(d)
	recode keep (0=1) if maxk<(121-obs) & maxk!=. & t>1961
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1966 & keep!=1
	egen k=group(tid) if t==1966
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1966
	egen obs=sum(d)
	recode keep (0=1) if maxk<(136-obs) & maxk!=. & t>1965
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1971 & keep!=1
	egen k=group(tid) if t==1971
	by id: egen maxk=max(k) 
	g d=1 if keep==1 & t==1971
	egen obs=sum(d)
	recode keep (0=1) if maxk<(141-obs) & maxk!=. & t>1970
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1975 & keep!=1
	egen k=group(tid) if t==1975
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1975
	egen obs=sum(d)
	recode keep (0=1) if maxk<(151-obs) & maxk!=. & t>1974
	drop k maxk tid d obs
	g tid=rnormal(0,100) if t==1991 & keep!=1
	egen k=group(tid) if t==1991
	
	by id: egen maxk=max(k)
	g d=1 if keep==1 & t==1991
	egen obs=sum(d)
	recode keep (0=1) if maxk<(171-obs) & maxk!=. & t>1990
	drop if keep!=1
	keep id t x lng 
	save "tempfinal.dta",replace
end




*Set up programs to simulate regression results based on simulated datasets*
//simulate results for main results, restricted time periods and number of countries, and specifications that include lagged dependent variable
program define psim,rclass
	local C = `1'
	//simulate panel dataset
	pdata
	xtset id t
	 g samp0=1 if t>1899 & t!=. 
	 g samp=1 if t>1969 & t!=. 
	 //run regular regression
    reghdfe lng  x , absorb(t id ) cluster(id)
	return scalar b1 = _b[x]
	return scalar s1 = _se[x]
	//restrict time period
	 reghdfe lng  x  if samp0==1,absorb(t id )  cluster(id)
	return scalar b2 = _b[x]
	return scalar s2 = _se[x]
	//restrict time period
	 reghdfe lng  x  if samp==1,absorb(t id )  cluster(id)
	return scalar b3 = _b[x]
	return scalar s3 = _se[x]
	//Include lagged dependent variable
  reghdfe lng  x l.lng , absorb(t id ) cluster(id)
	return scalar bg1 = _b[x]
	return scalar sg1 = _se[x]
	//restrict time period
	 reghdfe lng  x l.lng if samp0==1,absorb(t id )  cluster(id)
	return scalar bg2 = _b[x]
	return scalar sg2 = _se[x]
	//restrict time period
	 reghdfe lng  x l.lng  if samp==1,absorb(t id )  cluster(id)
	return scalar bg3 = _b[x]
	return scalar sg3 = _se[x]
	//drop random countries to restrict the sample size to 125 countries
bys id: gen idtag=_n==1
	bys idtag (t): gen rn=runiform()
	bys idtag (rn): gen select=(idtag==1) & (_n <= 125)
	bys id (t): replace select=sum(select)
	bys id (t): replace select=select[_N]
	keep if select==1
	  reghdfe lng  x , absorb(t id ) cluster(id)
	return scalar bc1 = _b[x]
	return scalar sc1 = _se[x]
	reghdfe lng  x l.lng, absorb(t id ) cluster(id)
	return scalar bgc1 = _b[x]
	return scalar sgc1 = _se[x]
	drop select idtag rn
		//drop random countries to restrict the sample size to 75 countries
	bys id: gen idtag=_n==1
	bys idtag (t): gen rn=runiform()
	bys idtag (rn): gen select=(idtag==1) & (_n <= 75)
	bys id (t): replace select=sum(select)
	bys id (t): replace select=select[_N]
	keep if select==1
	  reghdfe lng  x , absorb(t id ) cluster(id)
	return scalar bcc1 = _b[x]
	return scalar scc1 = _se[x]
	reghdfe lng  x l.lng, absorb(t id ) cluster(id)
	return scalar bgcc1 = _b[x]
	return scalar sgcc1 = _se[x]
end
//simulate results for interactions
program define psimint,rclass
	local C = `1'
	pdata
	sort id t
	g d=1
	//Create indicator that split countries into two groups and construct interaction variable
	by id: gen sum=sum(d)
		g zr=uniform() if sum==1
	egen zz=cut(zr),group(2)
	by id: egen z=max(zz)
	xtset id t
	g interaction=z*x
	g y=lng
	g y2=lng
	g y3=lng
	g y4=lng
	g y5=lng
	g y6=lng
	//Simulate different size interaction effects based on two main effects (0.098 and 0.16)
 replace y=lng+2*0.16*x if z==1
	  replace y2=lng+1.5*0.16*x if z==1
	   replace y2=lng+0.5*0.16*x if z==0
	   	  replace y3=lng+1.25*0.16*x if z==1
	   replace y3=lng+0.75*0.16*x if z==0
	    replace y4=lng+2*0.098*x if z==1
	  replace y5=lng+1.5*0.098*x if z==1
	   replace y5=lng+0.5*0.098*x if z==0
	   	  replace y6=lng+1.25*0.098*x if z==1
	   replace y6=lng+0.75*0.098*x if z==0
    reghdfe y x z interaction, absorb(t id) cluster(id)
	return scalar b1 = _b[interaction]
	return scalar s1 = _se[interaction]
	 reghdfe y2 x z interaction, absorb(t id) cluster(id)
	return scalar b2 = _b[interaction]
	return scalar s2 = _se[interaction]
	 reghdfe y3 x z interaction, absorb(t id) cluster(id)
	return scalar b3 = _b[interaction]
	return scalar s3 = _se[interaction]
	    reghdfe y4 x z interaction, absorb(t id) cluster(id)
	return scalar b4 = _b[interaction]
	return scalar s4 = _se[interaction]
	 reghdfe y5 x z interaction, absorb(t id) cluster(id)
	return scalar b5 = _b[interaction]
	return scalar s5 = _se[interaction]
	 reghdfe y6 x z interaction, absorb(t id) cluster(id)
	return scalar b6 = _b[interaction]
	return scalar s6 = _se[interaction]
end
//simulate results for dynamic effects
program define pdyn, rclass
	local C = `1'
	pdata
	xtset id t
		//Simulate different size dynamic effects 
   	g dyne=time_dem*0.0105
	replace dyne=0.21 if time_dem>20
	 	g dyne2=time_dem*0.0065
	replace dyne2=0.1303 if time_dem>20
		 	g dyne3=time_dem*0.00804
	replace dyne3=0.1609 if time_dem>20
	 g y=lng+dyne
	  g y2=lng+dyne2
	    g y3=lng+dyne3
	 tab time_dem,g(td)
	 tab time_pre,g(pre)
	 drop td1 pre1 pre2
    reghdfe y pre* td* , absorb(t id) cluster(id)
	return scalar b1=_b[td2] 
	return scalar s1=_se[td2]
		return scalar b2=_b[td3] 
	return scalar s2=_se[td3]
		return scalar b3=_b[td4] 
	return scalar s3=_se[td4]
		return scalar b4=_b[td5] 
	return scalar s4=_se[td5]
		return scalar b5=_b[td6] 
	return scalar s5=_se[td6]
		return scalar b6=_b[td7] 
	return scalar s6=_se[td7]
		return scalar b7=_b[td8] 
	return scalar s7=_se[td8]
		return scalar b8=_b[td9] 
	return scalar s8=_se[td9]
		return scalar b9=_b[td10] 
	return scalar s9=_se[td10]
			return scalar b10=_b[td11] 
	return scalar s10=_se[td11]
			return scalar b11=_b[td12] 
	return scalar s11=_se[td12]
			return scalar b12=_b[td13] 
	return scalar s12=_se[td13]
			return scalar b13=_b[td14] 
	return scalar s13=_se[td14]
			return scalar b14=_b[td15] 
	return scalar s14=_se[td15]
			return scalar b15=_b[td16] 
	return scalar s15=_se[td16]
			return scalar b16=_b[td17] 
	return scalar s16=_se[td17]
			return scalar b17=_b[td18] 
	return scalar s17=_se[td18]
			return scalar b18=_b[td19] 
	return scalar s18=_se[td19]
			return scalar b19=_b[td20] 
	return scalar s19=_se[td20]
			return scalar b20=_b[td21] 
	return scalar s20=_se[td21]
	
	    reghdfe y2 pre* td* , absorb(t id) cluster(id)
	return scalar bb1=_b[td2] 
	return scalar ss1=_se[td2]
		return scalar bb2=_b[td3] 
	return scalar ss2=_se[td3]
		return scalar bb3=_b[td4] 
	return scalar ss3=_se[td4]
		return scalar bb4=_b[td5] 
	return scalar ss4=_se[td5]
		return scalar bb5=_b[td6] 
	return scalar ss5=_se[td6]
		return scalar bb6=_b[td7] 
	return scalar ss6=_se[td7]
		return scalar bb7=_b[td8] 
	return scalar ss7=_se[td8]
		return scalar bb8=_b[td9] 
	return scalar ss8=_se[td9]
		return scalar bb9=_b[td10] 
	return scalar ss9=_se[td10]
			return scalar bb10=_b[td11] 
	return scalar ss10=_se[td11]
			return scalar bb11=_b[td12] 
	return scalar ss11=_se[td12]
			return scalar bb12=_b[td13] 
	return scalar ss12=_se[td13]
			return scalar bb13=_b[td14] 
	return scalar ss13=_se[td14]
			return scalar bb14=_b[td15] 
	return scalar ss14=_se[td15]
			return scalar bb15=_b[td16] 
	return scalar ss15=_se[td16]
			return scalar bb16=_b[td17] 
	return scalar ss16=_se[td17]
			return scalar bb17=_b[td18] 
	return scalar ss17=_se[td18]
			return scalar bb18=_b[td19] 
	return scalar ss18=_se[td19]
			return scalar bb19=_b[td20] 
	return scalar ss19=_se[td20]
			return scalar bb20=_b[td21] 
	return scalar ss20=_se[td21]
	
		    reghdfe y3 pre* td* , absorb(t id) cluster(id)
	return scalar bbb1=_b[td2] 
	return scalar sss1=_se[td2]
		return scalar bbb2=_b[td3] 
	return scalar sss2=_se[td3]
		return scalar bbb3=_b[td4] 
	return scalar sss3=_se[td4]
		return scalar bbb4=_b[td5] 
	return scalar sss4=_se[td5]
		return scalar bbb5=_b[td6] 
	return scalar sss5=_se[td6]
		return scalar bbb6=_b[td7] 
	return scalar sss6=_se[td7]
		return scalar bbb7=_b[td8] 
	return scalar sss7=_se[td8]
		return scalar bbb8=_b[td9] 
	return scalar sss8=_se[td9]
		return scalar bbb9=_b[td10] 
	return scalar sss9=_se[td10]
			return scalar bbb10=_b[td11] 
	return scalar sss10=_se[td11]
			return scalar bbb11=_b[td12] 
	return scalar sss11=_se[td12]
			return scalar bbb12=_b[td13] 
	return scalar sss12=_se[td13]
			return scalar bbb13=_b[td14] 
	return scalar sss13=_se[td14]
			return scalar bbb14=_b[td15] 
	return scalar sss14=_se[td15]
			return scalar bbb15=_b[td16] 
	return scalar sss15=_se[td16]
			return scalar bbb16=_b[td17] 
	return scalar sss16=_se[td17]
			return scalar bbb17=_b[td18] 
	return scalar sss17=_se[td18]
			return scalar bbb18=_b[td19] 
	return scalar sss18=_se[td19]
			return scalar bbb19=_b[td20] 
	return scalar sss19=_se[td20]
			return scalar bbb20=_b[td21] 
	return scalar sss20=_se[td21]
	
end
//simulate results for civil war onset as outcome
program define psimonset,rclass
	local C = `1'
	pdata2	
    reghdfe onset  x , absorb(t id ) cluster(id)
	return scalar b1 = _b[x]
	return scalar s1 = _se[x]
end
//simulate results with polyarchy index as treatment
program define ppoly,rclass
	local C = `1'
	pdata3
	xtset id t

    reghdfe lng  x , absorb(t id ) cluster(id)
	return scalar b1 = _b[x]
	return scalar s1 = _se[x]
end


*Run simulations*
simulate reg1_b=r(b1) reg1_s=r(s1) reg2_b=r(b2) reg2_s=r(s2) reg3_b=r(b3) reg3_s=r(s3) regg1_b=r(bg1) regg1_s=r(sg1) regg2_b=r(bg2) regg2_s=r(sg2) regg3_b=r(bg3) regg3_s=r(sg3) regc1_b=r(bc1) regc1_s=r(sc1) regcg1_b=r(bgc1) regcg1_s=r(sgc1) regcc1_b=r(bcc1) regcc1_s=r(scc1) regccg1_b=r(bgcc1) regccg1_s=r(sgcc1), reps(10000): psim 185
save "sim_0_time.dta",replace 

simulate reg1_b=r(b1) reg1_s=r(s1) reg2_b=r(b2) reg2_s=r(s2) reg3_b=r(b3) reg3_s=r(s3) reg4_b=r(b4) reg4_s=r(s4) reg5_b=r(b5) reg5_s=r(s5) reg6_b=r(b6) reg6_s=r(s6), reps(10000): psimint 185
save "sim_int.dta",replace 

simulate reg1_b=r(b1) reg1_s=r(s1), reps(10000): psimonset 185
save "sim_cw.dta",replace 

simulate reg1_b=r(b1) reg1_s=r(s1), reps(10000): ppoly 185
save "sim_poly.dta",replace 

simulate b1=r(b1) s1=r(s1) b2=r(b2) s2=r(s2) b3=r(b3) s3=r(s3) b4=r(b4) s4=r(s4) b5=r(b5) s5=r(s5) b6=r(b6) s6=r(s6) b7=r(b7) s7=r(s7) b8=r(b8) s8=r(s8) b9=r(b9) s9=r(s9) b10=r(b10) s10=r(s10) b11=r(b11) s11=r(s11) b12=r(b12) s12=r(s12) b13=r(b13) s13=r(s13) b14=r(b14) s14=r(s14) b15=r(b15) s15=r(s15) b16=r(b16) s16=r(s16) b17=r(b17) s17=r(s17) b18=r(b18) s18=r(s18) b19=r(b19) s19=r(s19) b20=r(b20) s20=r(s20) bb1=r(bb1) ss1=r(ss1) bb2=r(bb2) ss2=r(ss2) bb3=r(bb3) ss3=r(ss3) bb4=r(bb4) ss4=r(ss4) bb5=r(bb5) ss5=r(ss5) bb6=r(bb6) ss6=r(ss6) bb7=r(bb7) ss7=r(ss7) bb8=r(bb8) ss8=r(ss8) bb9=r(bb9) ss9=r(ss9) bb10=r(bb10) ss10=r(ss10) bb11=r(bb11) ss11=r(ss11) bb12=r(bb12) ss12=r(ss12) bb13=r(bb13) ss13=r(ss13) bb14=r(bb14) ss14=r(ss14) bb15=r(bb15) ss15=r(ss15) bb16=r(bb16) ss16=r(ss16) bb17=r(bb17) ss17=r(ss17) bb18=r(bb18) ss18=r(ss18) bb19=r(bb19) ss19=r(ss19) bb20=r(bb20) ss20=r(ss20) bbb1=r(bbb1) sss1=r(sss1) bbb2=r(bbb2) sss2=r(sss2) bbb3=r(bbb3) sss3=r(sss3) bbb4=r(bbb4) sss4=r(sss4) bbb5=r(bbb5) sss5=r(sss5) bbb6=r(bbb6) sss6=r(sss6) bbb7=r(bbb7) sss7=r(sss7) bbb8=r(bbb8) sss8=r(sss8) bbb9=r(bbb9) sss9=r(sss9) bbb10=r(bbb10) sss10=r(sss10) bbb11=r(bbb11) sss11=r(sss11) bbb12=r(bbb12) sss12=r(sss12) bbb13=r(bbb13) sss13=r(sss13) bbb14=r(bbb14) sss14=r(sss14) bbb15=r(bbb15) sss15=r(sss15) bbb16=r(bbb16) sss16=r(sss16) bbb17=r(bbb17) sss17=r(sss17) bbb18=r(bbb18) sss18=r(sss18) bbb19=r(bbb19) sss19=r(sss19) bbb20=r(bbb20) sss20=r(sss20) , reps(1000): pdyn 185
save "sim_dyn.dta",replace

***Produce figures***
*********************
use"sim_0_time.dta",replace
//calculate how the share of significant coefficients change as a function of effect size
g id=_n
expand 300
g d=0.001
bys id: gen sum=sum(d)
replace sum=0 if sum>0.29999
drop d
rename sum esize
replace reg1_b=reg1_b+esize
g tval=reg1_b/reg1_s
g sig=0
recode sig (0=1) if tval>1.96
bys esize: egen mean=mean(sig)
replace regg1_b=regg1_b+esize
g tval2=regg1_b/regg1_s
g sig2=0
recode sig2 (0=1) if tval2>1.96
bys esize: egen mean2=mean(sig2)
replace regc1_b=regc1_b+esize
g tvalc=regc1_b/regc1_s
g sigc=0
recode sigc (0=1) if tvalc>1.96
bys esize: egen meanc=mean(sigc)
replace regcg1_b=regcg1_b+esize
g tvalcg=regcg1_b/regcg1_s
g sigcg=0
recode sigcg (0=1) if tvalcg>1.96
bys esize: egen meancg=mean(sigcg)
replace regcc1_b=regcc1_b+esize
g tvalcc=regcc1_b/regcc1_s
g sigcc=0
recode sigcc (0=1) if tvalcc>1.96
bys esize: egen meancc=mean(sigcc)
replace regccg1_b=regccg1_b+esize
g tvalccg=regccg1_b/regccg1_s
g sigccg=0
recode sigccg (0=1) if tvalccg>1.96
bys esize: egen meanccg=mean(sigccg)
replace reg2_b=reg2_b+esize
g tvalt=reg2_b/reg2_s
g sigt=0
recode sigt (0=1) if tvalt>1.96
bys esize: egen meant=mean(sigt)
replace reg3_b=reg3_b+esize
g tvaltt=reg3_b/reg3_s
g sigtt=0
recode sigtt (0=1) if tvaltt>1.96
bys esize: egen meantt=mean(sigtt)
replace regg2_b=regg2_b+esize
g tvalgt=regg2_b/regg2_s
g siggt=0
recode siggt (0=1) if tvalgt>1.96
bys esize: egen meangt=mean(siggt)
replace regg3_b=regg3_b+esize
g tvalgtt=regg3_b/regg3_s
g siggtt=0
recode siggtt (0=1) if tvalgtt>1.96
bys esize: egen meangtt=mean(siggtt)

g gesize=esize*100

*Fig. 1: Main results*
capture graph drop *
tw (scatteri 1 0.051 1 0.132,recast(area) color(gray%20) lwidth(none)) (scatteri 1 0.088 1 0.09,recast(area) color(gray%60) lwidth(none)) (scatteri 1 0.12 1 0.122,recast(area) color(gray) lwidth(none)) (scatteri 1 0.149 1 0.151,recast(area) color(black) lwidth(none))   (line mean esize,lcolor(black))   if esize<0.251,yline(0.8,lcolor(black) lpattern(dash)) yline(0.9,lcolor(black)) legend(order(1 "IQR MV" 2 "Avg. MV" 3 "CRM 2020"  4  "ANRR 2019" 5 "Share sig.") ) graphregion(color(white)) title("Level",color(black)) xtitle(Effect size) ytitle(Share significant) name(a)
tw (scatteri 1 0.29 1 0.57,recast(area) color(gray%20) lwidth(none)) (scatteri 1 0.413 1 0.425,recast(area) color(gray%60) lwidth(none))  (scatteri 1 0.435 1 0.447,recast(area) color(black) lwidth(none))   (line mean2 gesize,lcolor(black))   if gesize<1.1,yline(0.8,lcolor(black) lpattern(dash)) yline(0.9,lcolor(black)) legend(order(1 "IQR MV" 2 "Avg. MV" 3 "KW 2015" 4 "Share sig.") ) graphregion(color(white)) title("Growth (in pct)",color(black)) xtitle(Effect size) ytitle(Share significant) name(b)
graph combine a b, graphregion(color(white))
graph export "main.pdf",replace

*Fig. 2: Limiting the number of countries*
capture graph drop *
tw (scatteri 1 0.051 1 0.132,recast(area) color(gray%20) lwidth(none)) (scatteri 1 0.088 1 0.09,recast(area) color(gray%60) lwidth(none)) (scatteri 1 0.12 1 0.122,recast(area) color(gray) lwidth(none)) (scatteri 1 0.149 1 0.151,recast(area) color(black) lwidth(none))   (line mean esize,lcolor(black)) (line meanc esize,lcolor(black) lpattern(dash)) (line meancc esize,lcolor(gray))  if esize<0.251,yline(0.8,lcolor(black) lpattern(dash)) yline(0.9,lcolor(black)) legend(order(1 "IQR MV" 2 "Avg. MV" 3 "CRM 2020"  4  "ANRR 2019" 5 "Full sample" 6 "125 countries" 7 "75 countries") ) graphregion(color(white)) title("Level",color(black)) xtitle(Effect size) ytitle(Share significant) name(a)
tw (scatteri 1 0.29 1 0.57,recast(area) color(gray%20) lwidth(none)) (scatteri 1 0.413 1 0.425,recast(area) color(gray%60) lwidth(none))  (scatteri 1 0.435 1 0.447,recast(area) color(black) lwidth(none))   (line mean2 gesize,lcolor(black)) (line meancg gesize,lcolor(black) lpattern(dash)) (line meanccg gesize,lcolor(gray))  if gesize<1.1,yline(0.8,lcolor(black) lpattern(dash)) yline(0.9,lcolor(black)) legend(order(1 "IQR MV" 2 "Avg. MV" 3 "KW 2015" 4 "Full sample" 5 "125 countries" 6 "75 countries") ) graphregion(color(white)) title("Growth (in pct)",color(black)) xtitle(Effect size) ytitle(Share significant) name(b)
graph combine a b, graphregion(color(white))
graph export "climit.pdf",replace

*Fig. 4: Limiting time periods*
capture graph drop *
tw (scatteri 1 0.051 1 0.132,recast(area) color(gray%20) lwidth(none)) (scatteri 1 0.088 1 0.09,recast(area) color(gray%60) lwidth(none)) (scatteri 1 0.12 1 0.122,recast(area) color(gray) lwidth(none)) (scatteri 1 0.149 1 0.151,recast(area) color(black) lwidth(none))   (line mean esize,lcolor(black)) (line meant esize,lcolor(black) lpattern(dash)) (line meantt esize,lcolor(gray))  if esize<0.251,yline(0.8,lcolor(black) lpattern(dash)) yline(0.9,lcolor(black)) legend(order(1 "IQR MV" 2 "Avg. MV" 3 "CRM 2020"  4  "ANRR 2019" 5 "Full sample" 6 "1900-" 7 "1970-") ) graphregion(color(white)) title("Level",color(black)) xtitle(Effect size) ytitle(Share significant) name(a)
tw (scatteri 1 0.29 1 0.57,recast(area) color(gray%20) lwidth(none)) (scatteri 1 0.413 1 0.425,recast(area) color(gray%60) lwidth(none))  (scatteri 1 0.435 1 0.447,recast(area) color(black) lwidth(none))   (line mean2 gesize,lcolor(black)) (line meancg gesize,lcolor(black) lpattern(dash)) (line meanccg gesize,lcolor(gray))  if gesize<1.1,yline(0.8,lcolor(black) lpattern(dash)) yline(0.9,lcolor(black)) legend(order(1 "IQR MV" 2 "Avg. MV" 3 "KW 2015" 4 "Full sample" 5 "1900-" 6 "1970-") ) graphregion(color(white)) title("Growth (in pct)",color(black)) xtitle(Effect size) ytitle(Share significant) name(b)
graph combine a b, graphregion(color(white))
graph export "tlimit.pdf",replace

*Fig. 3: Time periods included and transition events in sample*
capture graph drop *
use "democracy_for_simu.dta" ,replace
rename id id2
merge 1:1 id2 t using "gdp_for_simu.dta"
keep if _merge==3
drop _merge

xtset id t
g ch=abs(x-l.x)
sort id t
by id: gen ch2=sum(ch)
g ch3=0
recode ch3 (0=1) if ch2==1 & ch==1
collapse (sum) ch ch3,by(t)
g revt=t*-1
sort revt
gen cumevents=sum(ch)
gen fevents=sum(ch3)
g cum2=cum/225
twoway (bar fevents t if t>1849, fcolor(gray) lwidth(none)), graphregion(color(white)) xlabel(1850 (20) 2020) xtitle(Sample cut-off year,size(small)) ytitle(Number of countries with transition event in sample,size(small)) name(aa)
twoway (area cum2 t if t>1849, fcolor(bluishgray) lwidth(none)), graphregion(color(white)) xlabel(1850 (20) 2020) xtitle(Sample cut-off year,size(small)) ytitle(Share of transition events included in sample,size(small)) name(a)

use "democracy_for_simu.dta" ,replace
rename id id2
merge 1:1 id2 t using "gdp_for_simu.dta"
keep if _merge==3
drop _merge
rename id2 id
rename t year
merge 1:1 id year using "id_w_continent.dta"
keep if _merge==3
drop _merge
rename year t
xtset id t
g ch=abs(x-l.x)

collapse (sum) ch,by(t contine)
g revt=t*-1


sort cont revt
by cont: gen cumevents=sum(ch)
bys t: egen sum=sum(cumevents)

replace cum=cum/sum
twoway (line cum t if t>1849 & cont==0, lcolor(black)) (line cum t if t>1849 & cont==1, lcolor(black) lpattern(dash)) (line cum t if t>1849 & cont==2,lcolor(black) lpattern(dot)) (line cum t if t>1849 & cont==3,lcolor(gray)) (line cum t if t>1849 & cont==4,lcolor(gray) lpattern(dash))  if t<2011, graphregion(color(white)) xlabel(1850 (20) 2020) xtitle(Sample cut-off year,size(small)) ytitle(Share of evets in sample located in each continent,size(small)) name(b)  legend(order(1 "Europe" 2 "Africa" 3 "Asia" 4 "Oceania/Caribbean" 5 "Americas") region(lstyle(none) lcolor(white)) size(*1) symxsize(*.4) symysize(*.4) forcesize ring(0) position(11) colgap(tiny) row(1)) ylabel(0 (0.2) 0.6)


graph combine aa b,graphregion(color(white)) rows(2) cols(1)
graph export "slug.pdf",replace


*Fig. 5: Dynamic effects*
use "sim_dyn.dta",replace 
g id=_n
expand 20
g d=1
bys id: gen sum=sum(d)
drop d
foreach i in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 {
g t`i'=b`i'/s`i'
g tt`i'=bb`i'/ss`i'
g ttt`i'=bbb`i'/sss`i'
capture g s=0 
recode s (0=1) if t`i'>1.96 & t`i'!=. & sum==`i'
capture g sss=0 
recode sss (0=1) if ttt`i'>1.96 & ttt`i'!=. & sum==`i'
capture g ss=0 
recode ss (0=1) if tt`i'>1.96 & tt`i'!=. & sum==`i'
}

collapse (mean) s ss sss,by(sum)
capture graph drop *
tw  (line s sum,lcolor(black)) (line sss sum,lcolor(black) lpattern(dash)) (line ss sum,lcolor(gray)) ,yline(0.8,lcolor(black) lpattern(dash)) yline(0.9,lcolor(black)) legend(order(1 "MDE" 2 "CRM 2020" 3 "MV") ) graphregion(color(white)) title("") xtitle("Years since democratization") ytitle(Share significant)  ytitle(Share significant) ylabel(0 (0.2) 1) xlabel(1 (2) 20)
graph export "dynamic.pdf",replace


*Fig. 6: Interactions*
use "sim_int.dta",replace 
g id=_n
expand 3
g d=1
bys id: gen sum=sum(d)
drop d
g tv1=reg1_b/reg1_s
g sig1=0
recode sig1 (0=1) if tv1>1.96
g tv2=reg2_b/reg2_s
g sig2=0
recode sig2 (0=1) if tv2>1.96
g tv3=reg3_b/reg3_s
g sig3=0
recode sig3 (0=1) if tv3>1.96
g tv4=reg4_b/reg4_s
g sig4=0
recode sig4 (0=1) if tv4>1.96
g tv5=reg5_b/reg5_s
g sig5=0
recode sig5 (0=1) if tv5>1.96
g tv6=reg6_b/reg6_s
g sig6=0
recode sig6 (0=1) if tv6>1.96
g s1=.
g s2=.
replace s1=sig1 if sum==3
replace s1=sig2 if sum==2
replace s1=sig3 if sum==1
replace s2=sig4 if sum==3
replace s2=sig5 if sum==2
replace s2=sig6 if sum==1
collapse (mean) s1 s2,by(sum)
capture graph drop *
tw  (scatter s1 sum,mcolor(black)) (scatter s2 sum,mcolor(black) msymbol(circle_hollow))  ,yline(0.8,lcolor(black) lpattern(dash)) yline(0.9,lcolor(black)) legend(order(1 "Baseline effect=0.16 (MDE)" 2 "Baseline effect=0.098 (MV)") ) graphregion(color(white)) title("") xtitle("% difference in effect across groups") ytitle(Share significant) xlabel(1 "50%" 2 "100%" 3 "200%") ytitle(Share of interaction terms significant) ylabel(0 (0.2) 1)
graph export "interaction.pdf",replace


*Fig. 7: Consequences of low power*
use "sim_0_time.dta",replace
g id=_n
expand 300
g d=0.001
bys id: gen sum=sum(d)
replace sum=0 if sum==0.31
drop d
rename sum esize
replace reg1_b=reg1_b+esize
g tval=reg1_b/reg1_s
g sig=0
recode sig (0=1) if tval>1.96 | tval<-1.96
keep if sig==1
g wrong=0
recode wrong (0=1) if reg1_b<0
g exag=abs(reg1_b)/esize
bys esize: egen ws=mean(wrong)
bys esize: egen eg=mean(exag)
capture graph drop *
tw (line ws esize,lcolor(black)) if esize<0.2 & esize>0.02,graphregion(color(white)) name(a) title("Type S-Error",color(black)) ytitle(Share of estimates with wrong sign) xtitle(Effect size)
tw (line eg esize,lcolor(black)) if esize<0.2 & esize>0.02,graphregion(color(white)) name(b) title("Type M-Error",color(black)) ytitle(Exaggeration ratio) xtitle(Effect size)
graph combine b a,graphregion(color(white))
graph export "cons.pdf",replace

*Fig. A1: Polyarchy index as treatment*
use "sim_poly.dta",replace
g id=_n
expand 300
g d=0.001
bys id: gen sum=sum(d)
replace sum=0 if sum>0.29999
drop d
rename sum esize
replace reg1_b=reg1_b+esize
g tval=reg1_b/reg1_s
g sig=0
recode sig (0=1) if tval>1.96
bys esize: egen mean=mean(sig)
tw  (scatteri 1 0.058 1 0.068,recast(area) color(black) lwidth(none)) (scatteri 1 0.029 1 0.129,recast(area) color(gray%20) lwidth(none))  (line mean esize,lcolor(black))   if esize<0.251,yline(0.8,lcolor(black) lpattern(dash)) yline(0.9,lcolor(black)) legend(order(1 "MV avg." 2 "IQR MV" 3 "Share sig.") ) graphregion(color(white)) title("",color(black)) xtitle(Effect size) ytitle(Share significant) 
graph export "poly.pdf",replace

*Fig. A2: Civil war as outcome*
use "sim_cw.dta",replace
g id=_n
expand 300
g d=0.001
bys id: gen sum=sum(d)
replace sum=0 if sum>0.29999
drop d
rename sum esize
replace reg1_b=reg1_b+esize
g tval=reg1_b/reg1_s
g sig=0
recode sig (0=1) if tval>1.96
bys esize: egen mean=mean(sig)
tw  (scatteri 1 0.0098 1 0.0102,recast(area) color(black) lwidth(none))   (line mean esize,lcolor(black))   if esize<0.0251,yline(0.8,lcolor(black) lpattern(dash)) yline(0.9,lcolor(black)) legend(order(1 "BS 2018" 2 "Share sig.") ) graphregion(color(white)) title("",color(black)) xtitle(Effect size) ytitle(Share significant) 
graph export "cw.pdf",replace

***Run multiverse analyses***
*****************************
use "vpower.dta",replace
levelsof e_regiongeo, local(levels) 

 foreach i in lex_dem vdem_dem bmr_dem polity_dem fh_dem {

 foreach j in 1785 1899 1949 1969  {
 	foreach k in constant v2stfisccap e_pop {
	
  foreach l of local levels {
  	 use "vpower.dta",replace
 rename `i' democracy
xtset country_id year
 capture g constant=1
 reghdfe lng democracy `k'  if year>`j' & e_regiongeo!=`l',absorb(country_id year) cluster(country_id)	
  regsave democracy using "reg_`i'`k'`j'`l'.dta",level(95) ci addlabel(demindicator, `i',period, `j',control, `k',ex_region, `l') replace
 }
 }
 }
 }	                   
  	use "vpower.dta",replace
 	levelsof e_regiongeo, local(levels) 
	keep if _n==1
	                                      
 foreach i in lex_dem vdem_dem bmr_dem polity_dem fh_dem {

 foreach j in 1785 1899 1949 1969  {
 	foreach k in constant v2stfisccap e_pop {
	
  foreach l of local levels {
  capture append using "reg_`i'`k'`j'`l'.dta"
  capture erase "reg_`i'`k'`j'`l'.dta"
  }
	}
 }
 }
 save "multiverse_results.dta",replace
 //get average of estimates and interquartile range of estimates
 sum coef,detail
 
//run multiverse analysis with lagged dependent variable
use "vpower.dta",replace
levelsof e_regiongeo, local(levels) 

 foreach i in lex_dem vdem_dem bmr_dem polity_dem fh_dem {

 foreach j in 1785 1899 1949 1969  {
 	foreach k in constant v2stfisccap e_pop {
	
  foreach l of local levels {
  	 use "vpower.dta",replace
 rename `i' democracy
xtset country_id year
 capture g constant=1
 reghdfe lng democracy `k' l.lng if year>`j' & e_regiongeo!=`l',absorb(country_id year) cluster(country_id)	
  regsave democracy using "reg_`i'`k'`j'`l'.dta",level(95) ci addlabel(demindicator, `i',period, `j',control, `k',ex_region, `l') replace
 }
 }
 }
 }	                   
  	use "vpower.dta",replace
 	levelsof e_regiongeo, local(levels) 
	keep if _n==1
	                                      
 foreach i in lex_dem vdem_dem bmr_dem polity_dem fh_dem {

 foreach j in 1785 1899 1949 1969  {
 	foreach k in constant v2stfisccap e_pop {
	
  foreach l of local levels {
  capture append using "reg_`i'`k'`j'`l'.dta"
  capture erase "reg_`i'`k'`j'`l'.dta"
  }
	}
 }
 }
 save "multiverse_results_growth.dta",replace
 //get average of estimates and interquartile range of estimates
 sum coef,detail
 
//run multiverse analysis with polyarchy index as treatment
use "vpower.dta",replace
 	levelsof e_regiongeo, local(levels) 

 foreach i in poly {

 foreach j in 1785 1899 1949 1969  {
 	foreach k in constant v2stfisccap e_pop {
	
  foreach l of local levels {
  	 use "vpower.dta",replace
 rename `i' democracy

 capture g constant=1
 reghdfe lng democracy `k' if year>`j' & e_regiongeo!=`l',absorb(country_id year) cluster(country_id)	
  regsave democracy using "reg_`i'`k'`j'`l'.dta",level(95) ci addlabel(demindicator, `i',period, `j',control, `k',ex_region, `l') replace
 }
 }
 }
 }	                    
  	use "vpower.dta",replace
 	levelsof e_regiongeo, local(levels) 
	keep if _n==1
	                                      
 foreach i in poly {

 foreach j in 1785 1899 1949 1969  {
 	foreach k in constant v2stfisccap e_pop {
	
  foreach l of local levels {
  capture append using "reg_`i'`k'`j'`l'.dta"
  capture erase "reg_`i'`k'`j'`l'.dta"
  }
	}
 }
 }	
 save "multiverse_results_poly.dta",replace
 //get average of estimates and interquartile range of estimates
 sum coef,detail

